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Abstract 

Gravitational collapse of the cylindrical elongated cloud is studied by numerical magnetohydrody- 
namical simulations. In the infinitely long cloud in hydrostatic configuration, small perturbations 
grow by the gravitational instability. The most unstable mode indicated by a linear perturbation 
theory grows selectively even from a white noise. The growth rate agrees with that calculated by 
the linear theory. First, the density-enhanced region has an elongated shape, i.e., prolate spheroidal 
shape. As the collapse proceeds, the high-density fragment begins to contract mainly along the 
symmetry axis. Finally, a spherical core is formed in the non-magnetized cloud. In contrast, an 
oblate spheroidal dense disk is formed in a cloud in which the magnetic pressure is nearly equal 
to the thermal one. The radial size of the disk becomes proportional to the initial characteristic 
density scale-height in the r-direction. As the collapse proceeds, a slowly contracting dense part is 
formed (< 10% in mass) inside of the fast contracting disk. And this is separated from other part 
of the disk whose inflow velocity is accelerated as reaching the center of the core. From arguments 
on the Jeans mass and the magnetic critical mass, it is concluded that the fragments formed in 
a cylindrical elongated cloud can not be supported against the self-gravity and it will eventually 
collapse. 

Subject Headings: Interstellar: Matter — Interstellar: Magnetic Fields — Hydromagnetics - 
Stars: Formation 

1 Introduction 

The process of star formation begins in the interstellar clouds as a fragmentation of the clouds. 
Massive stars are formed only in giant molecular clouds, while less massive stars are born also in less 
massive dark clouds (for a review, see Larson 1991). This observational fact seems to suggest that the 
process of massive star formation and that of less massive stars are different. This seems to be related 
to the difference between the collapse of supercritical clouds and subcritical clouds (Shu, Adams, & 
Lizano 1987). If a mass of the cloud is larger than a critical mass (for a magnetic cloud, ~ magnetic 
flux/G 1 / 2 ), there is no equilibrium for the thermal pressure, the Lorentz force, and the centrifugal force 
to counter-balance the self-gravity. The supercritical clouds begin dynamical collapse. However, the 
subcritical magnetohydrostatic cloud is thought to evolve slowly only by the plasma drift (ambipolar 
diffusion; Mouschovias 1977, Nakano 1979) and/or the magnetic braking (Mouschovias 1979). The 
structure of the cloud changes in a relatively long time scale of the plasma drift and the magnetic 
braking. Finally, if the cloud becomes supercritical due to the decrease of magnetic flux at the center 
of the cloud or the decrease of angular momentum to support the cloud, it begins to dynamical 
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collapse, too. However, this is true only when the cloud in a static equilibrium is stable. When the 
time scale of some dynamical instabilities is shorter than the evolutionary time scale, the cloud may 
begin a dynamical contraction or fragmentation before it reaches the condition for the supercritical 
cloud. In a series papers, the fragmentation process in a subcritical, magnetohydrostatic cloud is 
studied by a non-linear magnetohydrodynamical simulation in two-dimension. This paper is devoted 
to the filamentary cloud and its gravitational fragmentation. 

The process of fragmentation and thus gravitational collapse of the high-density portion of the 
cloud are studied by a linear and non-linear analyses. Linear analyses of the gravitational instability 
in an isothermal slab have been done many authors (see Spitzer 1978 for a review). The qualita- 
tive answer is as follows: the slab is unstable for a perturbation whose wavelength is longer than a 
critical length and the most unstable perturbation has the wavelength of A max — 
and the e-growing time scale of T max ~ 2(47rG/o c ) _1//2 , where p c represents the density on the mid- 
plane of the disk (Elmegreen & Elmegreen 1978). The cylindrical cloud has a similar characteristic 
wavelength, A max — 20c s /(47rG / o c ) 1 ' /2 and a longer growth time scale as r max ~ 3(4vrGp c )- 1 / 2 (Na- 
gasawa 1987). It is concluded that if there is an inhomogeneity with an amplitude of Ap/p ~ 5%, 
in 3r max <10(47rGp c ) -1 / 2 the dense part of the cloud grows as a fragment with non-linear density 
contrast Ap/p ~ 1. This time scale is comparable to that for the ambipolar diffusion (Nakano 1988). 
Therefore, the dynamical instability in the subcritical static cloud should be studied. 

The isothermal non-magnetic cylindrical cloud has a critical mass per unit length beyond which 
no equilibrium is achieved as 
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where c s represents the isothermal sound speed in the cloud. The cylindrical cloud with a mass per 
unit length, A = A c , has the infinite density contrast between the center and the surface. As for the 
spherical isothermal cloud, the situation is different; the critical cloud has a finite density contrast of 
Pel Ps = 14 (Bonner 1956; Ebert 1955). The non-magnetic critical mass is equal to 

M c = 1.18 1/2 C * , (1.2) 

where p ex t represents the external pressure on the surface of the cloud. Finally, in the slab geometry, 
no critical column density of the self-gravitating disk exists. These differences come from the geometry 
or the dimension of the system. Generally, the magnetic field and the centrifugal force in rotating 
clouds have an effect of increasing the critical mass. From studying magnetohydrostatic equilibrium, 
the critical mass of the magnetized cloud is obtained as 
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where dm/d^B/G 1 / 2 ) means the mass-to-magnetic flux ratio at the center of the cloud (Tomisaka, 
Ikeuchi, &Nakamura 1988). Further the rotating cloud has a larger critical mass as 
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where M c mag and j represent, respectively, the mass which can be supported without any rotation 
[Eq. Ql.3] )] and a specific angular momentum of the cloud (Tomisaka, Ikeuchi, & Nakamura 1989). 
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Anyhow, the cloud with A > A c or M > M c eventually collapses, unless the excessive mass is eroded 
by any processes. 

Study of the dynamical evolution of the magnetized cloud has been restricted for supercritical 
clouds (Scott & Black 1980; Black & Scott 1982; Phillips 1986a, b; Dorfi 1982, 1989). The supercritical 
cloud contracts as a whole and forms a dense contracting core inevitably. The authors reported that 
they did not observed any fragmentation in the process of contraction. This seems to correspond to 
the fact that free-fall time scale is shorter than the growing time scale of the gravitational instability 
(see above). The situation seems to be much different for subcritical clouds. 

As for the non-magnetic cylindrical clouds, Bastien (1983), and Bastien et al. (1991) have studied 
the contraction of initially uniform cylindrical cloud with finite length. They achieved some conclu- 
sions on the fate of the above clouds: for example, for the elongated cloud with the ratio of length 
to diameter of the cylinder >2, using the initial Jeans number, Jo, which is defined as the ratio of 
gravitational to the thermal energies, the evolution is determined. When Jo^J2-frag^ t wo subconden- 
sations are formed. For more large Jeans number as J2-frag^>Jo^Jspindie, these two subcondensations 
collapses into one object. For the extreme case Jo^,J sp indiei the cloud contracts onto a line and forms 
a spindle. However, their simulation has a restriction that the initial state is far from the hydrostatic 
equilibrium. Excess free energy is liberated in the process of the relaxation from the initial state to 
the equilibrium, which may act an important role. To understand the relatively slow evolution in the 
subcritical cloud, we should take a (magneto-)hydrostatic configuration as the initial state for the 
simulation. 

2 Model and Numerical Method 
2.1 Initial Condition 

As described in the preceding section, we assume that the initial state is in a hydrostatic equilibrium. 
Using the cylindrical coordinate (z, r, (ft), if the initial cylindrical isothermal cloud is homogeneous 
in the z— direction, the magnetohydrostatic equilibrium configuration is obtained by the equation of 
hydrostatic balance and the Poisson equation as 



dift cjdp 1 dB< 
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ll(^)=4nG P , (2.2, 
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where ip, p, B, c s , and G represent, respectively, the gravitational potential, density, magnetic field, 
isothermal sound speed, and the gravitational constant. To derive equation( |2.1| ), it is assumed there 
is no helical magnetic field component {B^ = 0). The density distribution depends upon that of 
the magnetic field B z . In this paper, we restrict ourselves to two cases: (i) uniform magnetic field 
B z =constant, and (ii) the ratio of the magnetic pressure to the thermal one is constant. In the first 
model, the density distribution becomes identical to that of nonmagnetic isothermal cylinder as 

P{r) = ~ , 2 , (2-3) 

1 + ^^ 

c s 

where p c means the density at the center of the cylinder (r = 0). As easily seen, the density reaches 
zero only at the infinity (r = oo). Thus the isothermal cloud should be bounded by the external 
pressure, p ex t- The cloud has a boundary where p(r s ) is equal to p ex t/c^. The cloud radius is 
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expressed as 
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where p s is the density on the cloud surface and equals to p ex t/c 2 . In the present paper, we use the 
normalization as c s = AmG = p ey ± = 1. Thus the unit of the distance is chosen as H = c s / (A-nGp s ) 1 ^ 2 . 
The normalized density distribution, f(r) = p(r)/p a , and the radius of the surface, £ s = r s /H, are 
expressed as 

/(0 = (2.5) 



Fc2 




F l/2 _ ! 



1/2 



(2.6) 



where £ = r/H represents the normalized distance and F = p c /p s denotes the density contrast 
between the center and the surface. 

When the magnetic field plays a role in supporting the cloud, the density distribution is different 
from equation (p. Assuming that the magnetic pressure is proportional to the thermal pressure, 
i.e., B 2 /8ir oc c 2 p, the distribution of the density becomes as 
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where the parameter a is defined as 
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in terms of the plasma (3. This shows that the cloud becomes thick with increasing a, and the density 
scale-height changes in proportion to (1 + a/2) 1 / 2 . Typical density distributions are shown in Figure 
1. 

We assumed a hypothetical situation that the cloud is confined in a low-density ambient medium 
which has no importance as the source of gravity but has a finite pressure p ex t- Thus, the initial 
solution of p(r) has a boundary where p(r) = p cxt /c 2 . Beyond this cloud boundary, a tenuous and 
thus hot medium is assumed to extend. However, as seen in the next subsection, since we assume the 
isothermal equation of state, it is difficult to calculate two kinds of gases, cold cloud component and 
warm intercloud component, separately. Therefore, we adopt a one-fluid approximation, where the 
gas has an identical temperature; but the ambient gas is defined as the portion with p < p cxt /c 2 s = p s , 
and this gas has no effect of the gravitational field made by the cloud. This gives a virtual distribution 
of ambient gas as 

to counter-balance the gravity by the cloud as 



4(F 1 / 2 -l) 
F l/2 



(2.9) 



5(0 



dlnp 



(2.10) 



- 5 - 



2.2 Basic Equations 

The cylindrical symmetry is assumed: d/d(fi = 0. For the dense gas found in interstellar clouds with 
~ 10K, the equation of state is well approximated with the isothermal one. We assume here that 
the gas obeys the isothermal equation of state. Thus the basic equations become the equations of 
magnetohydrodynamics for isothermal gases as 

dp d , . 1 d , , 

dpvz d 1 d 2 dp dif) If dB r dB z \ 

-qT + Yz {pVzVz) + rd-r {rpVzVr) = - C *Tz ~ P ~dz ~ 4tt \dz d~r~) (2 ' 12) 

dpvr + d_, v s + lJL( rpv v \ = _ ( ?^P_ _ p ^t + — ( - 9Bz \ B (2 13) 

dt dz r dr s dr dr 4ir \ dz dr J ' 
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-JL = - — [r(v z B r - v r B z )], (2.14) 
dt r dr 

dB r d 



. }i —(v z B r -v r B z ), (2.15) 



dz 2 r dr \ dr 
where the variables have their ordinary meanings. 

To include only the gas density in the cloud as the source of the self-gravity, the right-hand side 
of the equation ( 2.16| ) is calculated as follows: 



P(z,r) = \ " 'i P>Ps > (2.17) 
1 ; \ 0, i£p<p t . K J 

As long as the cloud keeps its hydrostatic configuration, the solution of the above equations coincides 



with the initial configuration, equations (2.5) and ( |2.9| ). However, if the ambient matter accretes onto 
the cloud by the gravity or the cloud material evaporates by the pressure force, the above procedure 
may give wrong results. We monitored the mass of the cloud in numerical runs, and check the validity 
of the ambient matter condition. We confirmed that the cloud mass only changed less than ~ 3%. 

Nondimensional variables are used as follows: the density, p' = p/p s , the velocity, v' = v/c s , 
the gravitational potential, tp' = tp/c 2 , the magnetic fields, B' = B/(4-7rc^ / o s ) 1//2 , the time, t' = 
t(4irGp s ) l l 2 , and the linear size, r' = r/H. Using this normalization, equations ( 2.12[ ), ( p. 13 ), and 
(|2.2| ) become 

Other equations are identical even when these normalized variables are used. Hereafter, we will use 
these nondimensional variables and abbreviate the prime unless otherwise mentioned. 
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2.3 Numerical Method 

The basic equations are solved by the finite difference method. The mesh spacing (Az, Ar) is 
constant spatially except for models CB and CB2 in table 1. Unequal spacing is used for model CB 
and CB2 to see the fragment more closely. Cell numbers in one dimension are chosen N = 200 — 800 
(table 1). The "Monotonic Scheme" (van Leer 1977; Norman & Winkler 1986) is adopted to solve 
the hydrodynamical equations, i.e., equations ( |2.11| )-( ^.13| ), and the "Constrained Transport" method 
(Evans & Hawley 1988) is adopted to solve the induction equations of the magnetic field, i.e., equations 



(2.14) and (2.15). Although the equation of state is different, this scheme is essentially the same as 
that used by Tomisaka (1992), For numerical stability at the shock front, an artificial viscosity is 
included in equations ( |2.12j ) and ( |2.13 ). A tensor artificial viscosity term is used as shown in the 



Appendix of Tomisaka(1992). The program for adiabatic gas was checked by several test problems: 
Sedov solution, spherical stellar wind bubble solution (Weaver et al. 1977), one-dimensional magnetic 
Rieman problem (Brio & Wu 1988). Further, the contraction of a spherical isothermal cloud is 
compared with solutions of Larson (1969) and Penston (1969). The radial density distribution is well 
fitted by oc r~ 2 and derived velocity profiles agree with their results by one-dimensional calculation. 

To solve the Poisson equation is reduced to find the solution of simultaneous linear equations with 
a dimension of N 2 . The linear equations are solved by the "Conjugate Gradient Squared method" 
(CGS method, Dongarra et al. 1991) preconditioned by the modified incomplete LU decomposition 
(MILUCGS) (Meijerink & van der Vorst 1977; Gustafsson 1978). A brief description of this method 
and its performance are described in Appendix. 

On the upper (z = l z ) and lower [z = 0) boundary, the cyclic boundary condition is applied. On 
the outer boundary the fixed boundary condition is applied. 

3 Result 

3.1 A Typical Evolution 

First, the numerical result without any initial perturbation is seen. This corresponds to the model 
with a = 1 and F = c 2 p c /p ex t = p c /p s = 100 (model A). The size of the calculated region is taken as 
1.93H, which coincides with the wavelength which has the maximum growth rate in the gravitational 
instability (Nakamura, Hanawa, & Nakano 1993). Figure 2 shows the evolutions of density (left 
panel) and magnetic fields (right panel). In t = 3.6, no prominent fragmentation appears. In the 
stage shown in Figure 2c, a high-density region elongating in the z-direction appears. The shape 
of this high-density region is prolate, i.e., if it fits to an ellipsoid the major axis coincides with the 
symmetric axis (the z-axis). This agrees with the result of linear analyses (Nakamura et al. 1993; 
Nagasawa 1987); The wavelength of the most unstable perturbation was predicted from the linear 
analysis as ~ 2 x H. And the global shape of the high-density region in Figure 2c coincides with an 
expected shape from the eigen function of the linear analysis. From this figure, it is confirmed that 
the perturbation with a shorter wavelength than that of the most unstable mode does not appear 
at all, before the most unstable mode grows into the non-linear region. However, a perturbation 
with a longer wavelength than the numerical box can not be properly calculated. We have calculated 
in model A2 (not shown) the same problem with a larger spacing and thus with a numerical box 
twice as large as model A. In model A2, two high-density regions appear and each fragment shows 
the identical evolution as model A, which confirms the most unstable mode expected from the linear 
analysis grows first. 

After the stage in Figure 2c, collapse along the symmetry axis proceeds as well as the radial 
direction. The final structure is a disk contracting towards the center seen in Figure 2d. Since the 
gas easily falls down along the magnetic field lines, this disk is perpendicular to the symmetric axis 
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of the cylinder. The magnetic fields are squeezed by the effect of radial contraction near the disk 
(right panel). Since this configuration is unstable against the Parker instability (1979 and references 
therein), the gas is accumulated flowing along the magnetic fields. Contracting speed in the z- 
direction is ~ twice as fast as the radially contracting speed. In the central region, a high-density 
core is formed, on to which accretion of extra mass continues. Finally, near the high-density core 
the infalling speed exceeds the isothermal sound speed at this stage. It is concluded that a typical 
product of the self-gravitational instability in the cylindrical cloud is oblate spheroidal disks separated 
by A max - 

The phase of the position where the disk is formed (z c ) is not determined a priori. This is related 
to the question that where the irregularity comes from even if the initial state is uniform along the 
z-axis. After several numerical tests, it is concluded that the main contribution of irregularity is 
coming from the Poisson solver; Even if the density distribution, p, is uniform in the z-direction, 
the numerical values of the potential, tp, is not completely uniform. This is due to the CGS scheme. 
Although the typical relative error in the potential is small, Aip/ip ~ 1CP 7 , this irregularity grows by 
the effect of gravitational instability and finally it forms fragmentations. 

To confirm the idea, another model B is calculated, whose parameters are identical with model 
A except for the existence of initial irregularity. In model B, a perturbation with relative contrast 
5i n i = Ap/p = 10 -2 is added as 

/ 2ttz\ 
p(z,r) = p (r) l-<5 mi cos-— , (3.1) 



where po and l z represent, respectively, the initial density determined by /(£) of equation (|2.7| ) and 
the z-length of numerical box. The wavelength of this perturbation is again taken identical with that 
of the most unstable mode in the linear analysis. The structure at t = 1.26 is shown in Figure 3. This 
shows that the final structure of this model is identical with that of the model A (Fig. 2d), although 
the time necessary for the fragment to grow is shortened much (from t = 5.6 to t = 1.26). It is clear 
that this quick evolution is due to a finite amplitude of the perturbation added to the initial state, 
compared with the intrinsic irregularity from the Poisson solver in model A. However, the evolution 
itself is very similar with each other. 

To confirm this, the time evolution of the amplitude of the density irregularity is plotted for both 
models A and B in Figure 4. Using the average < p(r = 0) > and the maximum p(r = 0)| max density 
on the z-axis, the relative density contrast 

_ Pi r = 0)|max ., /„ „n 

< p[r = 0) > 

is plotted against the time passed after the calculation begins. This shows that in model A the 
perturbation grows in accordance with oc exp [ut] for l<i<5, from 10 -8 to 0.1. Model B shows the 
similar growth but the region in which it shows the exponential growth is much restricted. However, 
if the origin of time for model B is moved, two curves for models A and B coincides. This indicates 
that the evolution is similar whether the most unstable mode is overlaid on to the initial density 
distribution (model B) or not (model A). This confirms the expectation that the most unstable mode 
grows selectively even from the irregularity of white noise. Thus, we add the perturbation of the 
expected most unstable mode to the initial density distribution as model B. 

3.2 Non-Magnetic Cylinder 

To see the effect of the magnetic field, we calculated the non-magnetic isothermal cylinder as model C. 
From the linear analysis, the wavelength of the most unstable mode is equal to ~ 2.21H (Nakamura 
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et al. 1993). Since the magnetic field has an effect of supporting the cloud, the radial size of the cloud 
becomes more compact than that of models A and B. 

In Figure 5, the evolution of model C is shown. In t = 1.2, the fragment, whose density contrast 
5 reaches ~ 3, has entered in non-linear region (Fig. 5a). However, comparing with Figure 3 (these 
two models are assumed to have irregularities with the same initial amplitude), it is shown that this 
model evolves slower than model B. This difference is understood as follows: (1) since the magnetic 
fields of model B has a configuration unstable for the Parker instability, the system becomes more 
unstable including the magnetic fields, or (2) the mass per unit length of the cylinder and thus the 
effect of the self-gravity becomes larger when the magnetic fields are supporting the cloud. 

Figure 5b shows the structure at the age of t = 1.4. In this phase, the maximum density at the 
center of the core reaches 10 5 . The envelope with low-density ~ 10 1 ' 5 seems to have the shape of 
prolate spheroid as well as the magnetic models. In contrast, the high-density central region appears 
to be almost spherical symmetric. The difference between structures appearing in the magnetic cloud 
and the non- magnetic one is apparent. Since there is no lateral restoring force by the magnetic fields, 
the high density fragment shows spherical shape, while in the magnetic cloud the fragment is formed 
from the matter flowing along the magnetic fields as a disk perpendicular to the magnetic fields. 



3.3 Central Core 

To see the structure of the central core, the density profile is plotted against log \ z — z c \ and logr in 
Figure 6 for models B and C, where z c represents the position of the core center. 

There exist two characteristic power law solutions. One is a singular isothermal sphere solution 
(Chandrasekhar 1939): 



Ksmg 

or in a non-dimensional form: 
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The other is an asymptotic solution of equation fl2,5|) for £ S> 1 or r ^> H: 
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or in a non-dimensional form: 

Pout 1 



(3.6) 



As seen in Figure 6, the distribution in the r-direction is well fitted by equation (3.6) for the outer 
envelope. As collapse proceeds, the radial density distribution in the inner part of the cloud becomes 
similar to that of the singular isothermal sphere [equation (ft.3|)1 ■ 

As for the z-distribution, in the non-magnetic cloud (model C), the density distribution is very 
similar to that in the r-direction. However, in the magnetic cloud (model B), the r- and ^-distributions 
are much different. This difference comes from the effect of magnetic fields, that is, in models A and 
B a disk is finally formed, although in model C an almost spherical core is formed. Since the gravities 
by a disk and a sphere are quit different, the resultant density distributions are different. 

When the spherical cloud has the density distribution like equation (|3.3|) , the total mass inside of 
the sphere R is written as 

M(R) = j R p sing 4nr 2 dr = ^-R, 
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using the density p(R) at r = R. This represents the Jeans mass for the spherical cloud with infinite 
center-to-surface density contrast, Mj(R). Thus, when the density increases more rapid than equation 
( |3.3| ) toward the center, the actual mass is larger than this Jeans mass and contraction never stops. 
In contrast, when the density distribution is flat, the actual mass M(r) never surpasses Mj(r) in any 
region of r < R. Considering this, since the radial density distribution p^p s i n g, the spherical core 
formed in the center of the non-magnetic cloud seems to continue to contract. 

From the studies on magnetohydrostatic equilibrium (Tomisaka et al. 1988), the mass-to-flux 
ratio at the center of the cloud, M/^/G 1 / 2 ), has a crucial role to divide the clouds into super- and 
subcritical clouds; that is, considering a cloud whose mass is larger than the non-magnetic critical mass 
(J>ilf cr of eq.[ |1.2| ), when the mass-to-flux ratio is larger than ~ l/(2-7r) the cloud has no equilibrium 
solutions; while the cloud has at least one hydrostatic solution, when the ratio is smaller than this 
value. Since the mass-to-flux ratio remains constant along the contraction, the mass-to-flux ratio in 
the fragment can be written using quantities in the initial stage as 

M p c l z GV 2 pVXg^ . . 



Q/G 1 / 2 B c " (Ana) 1 / 



C.s 



where l z represents the z-length of the numerical box and is taken equal to the wavelength of the 
most unstable perturbation, ~ 20c s / '(4irG p c ) 1 ^ 2 ■ This ratio is approximately equal to 

M L59 (3.9) 



Thus, for the cloud with the mass appreciably larger than that of equation ( |l.2[) , only when the 
magnetic field is extremely strong as a>100, the magnetized fragment may be subcritical. The con- 
traction of the core can not be stopped by a magnetic field, which is assumed ordinarily, a ~ 1. This 
is valid until the wavelength of the most unstable perturbation is much longer than 20c s /(47rGp c ) 1 / 2 
in a cloud with low density contrast F. 

In conclusion, the cores formed in the fragment of the cylindrical cloud seems to continue to 
contract, until other effects not considered here may work in supporting the fragment, such as, the 
equation of state of the gas changes from the isothermal to more hard one. 

3.4 Diffuse Cloud 

Here, we compare the cloud with high central density and that with low density contrast calculated 
in models D. In Figure 7a, the evolution of model DB is plotted. Both models B and DB have the 
same parameters except for the center-to-surface density ratio, F. Difference between models B and 
DB are apparent: although a disk is formed for both models by the effect of the magnetic field, the 
disk perpendicular to the magnetic fields of model DB is much larger than that of model B. The disk 
has a radial size of 0.16 (model B) and 0.71 (model DB) (the size is measured as a radius of the 
largest contour line which has an oblate shape). 

The radial extent of the finally formed disk seems to correspond to the size of the initial central 
cylindrical core, i.e., in equation (|2.3j ) the initial density is almost uniform irrespective of the radial 
distance for r<r c = 2 L5 (1 + a) 1 / 2 c s /(47rGp c ) 1 / 2 . Since this gives 

r ~ 4 f^) ^ (1±*\ * (3 10) 



-10- 



the above result shows that a disk threaded by magnetic fields whose radial extent is equal to ~ 40 
- 55 % of the initial core size is formed. 

If this is true, does the cloud with a large a form a large disk? Yes. In model DD with a = 4 
(Fig. 7b), the disk with a radial size of ~ 1.2 is formed, while r c ~ 2(F/10) 1 / 2 [(1 + q)/5] 1 / 2 . Thus, it 
can be concluded that in the magnetized cylindrical cloud a disk is formed perpendicularly threaded 
by the magnetic fields and its radial extent is related to r c . 

The next characteristic point of the evolution of diffuse clouds is a low growth rate of fragmen- 
tation. A typical time scale for the perturbation to evolve into the nonlinear stage is t ~ 4.72 for 
a = (model DC), t ~ 4.12 for a = 1 (model DB), and t ~ 3.31 for a = 4 (model DD). This is much 
longer than that of models B and C with F = 100. This apparently comes from the fact that the 
growth time-scale is proportional to the free-fall time for the central density. Finally, comparing the 
non-magnetic clouds of models C and DC, Figures 4 (model C) and 7c (model DC) look very similar. 
These figures are scaled in proportional to l z , and l z is chosen as A max — 20c s /(47rGp c ) 1 / 2 . This 
means that the actual size of the spherical core of model DC is just 10 1 / 2 times larger than that of 
model C. If we scale the size by c s / (AnGpc) 1 / 2 instead of H = c s /(47rGp s ) 1 / 2 , the structures of these 
models becomes almost identical. Therefore, the structure of the fragment is mainly determined by 
the initial highest density at the center of the cylinder cloud for the non-magnetic cloud. 

3.5 Further Evolution 

To see the structure of the contacting central core in more detail, we have to execute larger calculation 
with finer resolution using a plenty of meshes. However, this approach will meet an inevitable limit 
of available memories on any computers. Therefore, we avoid this difficulty by varying mesh sizes, as 
a finer mesh is used around the expected position where the fragment is formed and a coarser mesh is 
used far from the cloud. We apply so-called log-mesh: the size of the grid is chosen like Azi = F z Azi-i 
for i > i c , Azi = F z Azi + \ for i < i c , and Arj = F r Arj-\ (F z , F r > 1), where i and j represent the 
sequential numbers of the cell in the z-direction and in the r-direction, respectively, and subscript c 
means the the midpoint of the z-axis. Models CB and CB2 have the same physical parameters as 
model B but they have ~ 2 - 4 times closer spatial resolution than model B, respectively. Although 
these two models use different initial amplitude of perturbations (table 1), qualitative evolutions are 
similar with each other. 

Due to a finer resolution, the evolution can be traced till the maximum density reaches ~ 10 6 = 
10 4 x initial central density. In Figure 8, we show the structure at t = 0.705 for model CB2 (due to 
a larger amplitude of <5i n ; = 10 _1 in model CB2, irregularities grow faster than models CB and B. 
This time scale is equivalent to t = 1.26 for models in which smaller initial perturbation is assumed 
as 5- m i = 10 -2 ). As shown in Figure 8c and d, collapse speed becomes faster as reaching the fragment 
center. And the maximum falling speed increases as collapse proceeds. From Figure 8c, at that time, 
the maximum infall velocity in the z-direction agrees with that observed in the asymptotic solution 
for the spherical isothermal collapse by Larson(1969) and Penston(1969), 3.28c s . It is shown that 
high-density region with a thickness of ~ 10 -2 , i.e., \z — z c \<10~ 2 H = 10 _1 c s /(47rGp c mit) 1 ^ 2 , is 
formed, in which contraction speed is decreased as reaching z c . In the r-direction, in contrast, very 
smooth distributions of density, magnetic field strength, and velocity are seen. This numerical result 
qualitatively agrees with the calculation by Hanawa, Nakamura, & Nakano (private communication) 
using a different numerical scheme. 

Thin dense object which is contracting slowly now seems to be separated from the outer inflow 
region seen in \z — z c |>10~ 2 . The fraction of mass which form this relatively static high-density disk 
is estimated as 6% - 10 %. So, if stars are formed from the relatively static matter, the ratio of the 
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mass of newly formed protostars to that of a parent cloud becomes as 6 % - 10 %. 

Therefore, the next evolution stage after this simulation ends seems to be a proto-stellar system 
consisting of a proto-star and an accreting matter onto the proto-star, whose luminosity comes from 
the liberated gravitational energy of falling material. 

4 Discussion 

4.1 A Model with Uniform Magnetic Fields 

In the preceding section, we investigated the models that the cloud is supported by magnetic fields. 
In contrast, here, we study the model with uniform magnetic field and compare the structure of 
fragments. Instead of B z /8ir oc c^p, we assume that the strength of B z is constant. Radial density 



distribution is given by equation (|2.3|) . Relative strength of the magnetic fields is expressed by a 
parameter a c which is equivalent to the ratio of B^/4ir to (? s p at the center of the cloud. 

In Figure 9, the structure at t = 1.74 is shown. This figure shows several differences with Figures 
1 and 2 (constant a model), (i) the cloud outer boundary: the distance from the z-axis increases near 
z ~ / z /2, and decreases far from the fragment. This seems to come from the fact that the magnetic 
fields becomes more important far from the z-axis and the gas moves toward z ~ l z /2 without any 
effects on the magnetic fields (Fig. 9b). As a result, the density contours indicates a convex shape. 
Contrarily, in models with constant a, the gas falls down along the valley of magnetic fields. Thus, 
the cylindrical tube of the cloud shrinks near z ~ l z /2 and inflates near z — and z ~ l z . (ii) the 
magnetic field line: as seen in Figure 9b the field line breaks weaker than previous models. This 
seems to come from the fact that the configuration is not unstable against the Parker instability, that 
is, the magnetic valley does not necessarily promote the instability in this case. That is, there is no 
positive feedback. The second factor is: since the ratio of magnetic pressure to the thermal one, a/2, 
increases outwardly, the magnetic field becomes relatively important for the dynamics with increasing 
r, even if a c is taken ~ 1. Thus the field line is almost locked in the initial configuration. 

We studied two models for the magnetic field configuration, a =constant and B z =constant, 
which correspond two extremes which are realized in the real interstellar space. 

4.2 Classification of the Cylindrical Cloud 

How does a cylindrical cloud with a finite line density A and a magnetic field strength evolve? In 
the preceding section, it is shown that even in a cloud in the hydrostatic equilibrium, fragments are 
formed and dense part of the fragment falls into a runaway collapse. From equation ( |2.7| ), the line 
density of a cloud with F = p c /ps, is written as 



2 C ? FV2 - 1 



a 



This leads to the maximum line density which can be supported by the thermal and magnetic pressure 
as 2 

A max = ^ffl + ^, for F^oo, (4.2) 



G V 2, 

which include the critical line density: equation( |l.l| ). 

Therefore, if the cloud has a line density smaller than A max , which is related to the isothermal 
sound speed, c s , and the Alfven speed, a 1//2 c s , the cloud has a magnetohydrostatic equilibrium. 
In such a cloud, irregularity grows in several xr max depending on the amplitude of irregularities. 
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Adopting 3r max as a typical growth time, this gives 

/ p n -1/2 / n -1/2 

r grow ~ 2.4 x 10 6 yr ^ t • (4.3) 

g y viooy \2 x 10- 22 g cm- 3 / V ' 

Final phase of the contraction, that is, the phase when the maximum density in the high-density core 
increases rapidly, continues only in a short duration ~0.1 x r grow . 

On the other hand, a cloud with super-critical line density A > A max does not have any magneto- 
hydrostatic configuration. Thus, the cloud contracts as a whole, and forms a thin cloud like a string. 
Recently, hydrodynamical simulation has been done for such a massive contracting cylindrical clouds 
by Inutsuka and Miyama (private communication) for non-magnetic clouds. They found that in such 
clouds, the irregularity grows relatively slowly compared with the contraction time-scale ~ free-fall 
time-scale and the cloud does not fragment but forms a string as long as the gas obeys the isothermal 
equation of state. 

Comparing these two, it is expected that super- and subcritical cylindrical cloud show completely 
different evolutions. This may be related to the mass function of the new-born stars. 

5 Summary 

We studied the process of fragmentation in an isothermal cylindrical cloud with infinite length. By 
a numerical magnetohydro dynamical method, the evolution of the fragment was investigated from 
the linear stage to nonlinear stage throughout. It is shown that the fastest growing perturbation has 
the same wavelength as predicted with the linear theory. In a linear stage the fragment appears as a 
prolate spheroidal shape. However, in a nonlinear stage, the fragment threaded by the magnetic field 
forms a disk perpendicular to the field line. At the center of the fragment there forms a high-density 
core, which continues to collapse. Non-magnetized cloud only forms collapsing spherical cores which 
are separated by A max . For magnetic clouds, it is shown that when the collapsing velocity reaches 
~ 3.5c s , relatively slowly contracting dense inner part is formed in the contracting disk, which also 
extends perpendicular to the magnetic fields. This disk is in almost static and is separated from the 
accretion flow which is accelerated as reaching the center of the core. 

I would like to thank T. Hanawa (Nagoya University) for stimulating discussion. We compared 
the numerical results with different schemes and could confirm the reliability of the schemes and 
the physics. I also thank S. Inutsuka, S. M. Miyama (National Astronomical Observatory) and 
M. Y. Fujimoto (Niigata University) for useful discussions. Numerical calculations were carried out 
by supercomputers: Hitac S-820/80's (Hokkaido University and University of Tokyo), Facom VP200 
(Institute of Space and Aeronautical Sciences). This work was supported in part by Grants-in-Aid 
for Science Research from the Ministry of Education, Science, and Culture (04233211 and 05217208). 
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Tables 



TABLE 1: Model Parameters 



Model 


a 


F 


<5ini 


h 


cell number 


A 


1 


100 





1.935 


400x400 


A2 


1 


100 





3.87 


400x400 


B 


1 


100 


10~ 2 


1.935 


400x400 


C 





100 


io- 2 


2.21 


400x400 


DB a ... 


1 


10 


io- 2 


1.935 


400x400 


DC a ... 





10 


io- 2 


2.21 


400x400 


DD a ... 


4 


10 


10~ 2 


1.935 


400x400 


CB 6 ... 


1 


100 


io- 2 


1.935 


400x280 


CB2 b .. 


1 


100 


io- 1 


1.935 


800x560 


UB C ... 


1 


100 


IO" 2 


1.935 


200x200 



a In a series of models with character D, diffuse clouds are considered. 

^Spatially varying spacing grid is used. The physical parameters are taken as identical as model B. 
c Uniform magnetic fields are assumed. In this model, a represents the value at the center of the cloud 
(r = 0), that is a c . 



TABLE Al: Necessary CPU time to solve the Poisson equation once 



CGS (sec) MILUCGS (sec) N 



Machine 



Sparc IPX (Workstation) 505 143 200 

M-682H (General Purpose Computer) 43.1 14.1 200 

S-820/80 (Supercomputer) 4.1 2A 400 
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A MILUCGS 



Here, a brief description of MILUCGS and its performance are described. The finite difference scheme 
for the Poisson equation is written 



/ 



ai 
db 2 



«2 



ubo 



UC 2 



\ ( Ik \ 

1p3 



dCm+l 



db m+ i a m+ i ub m+ i 



dc 



m+2 



V 



/ PI \ 

P2 
P3 



(A.l) 



\ Pnm ) 



where two-dimensional expression of tpij and pij (1 < i < m and 1 < j < n) are converted to a 
one-dimensional expression ^ and p/% with k = (j — 1) x m + i. m and n represent, respectively, 
the cell numbers in the z- and r-direction. The coefficients db and ub correspond to the derivatives 
with respect to the z-direction, and the coefficients of dc and uc correspond to those of r-direction. 
If the problem is a simple Neuman or Dirichret boundary condition, the matrix to be solved is like 



equation (A.l). However, since we apply here a cyclic boundary in the z-direction, we have to add two 
coefficients expressing this boundary condition. That is, the term of ip m j = ipjm should connect with 



m-l,j 



jm- 



-1 



«1 

db 2 



and ip 
ub\ 

«2 



hi 



ubo 



m—m+l 



ud^ 



instead of ipj rn+ i. The equation for cyclic boundary becomes 



ud 2 



\ 



uc 2 



ddm 

dc m+ \ dd m+ i db m+ i a m+ i ub m+ i 

dc m+2 dd m+2 



V 



■01 
1p2 



m+1 



Pi 

P2 



Pm+1 



\ Pnm J 



(A.2) 



where only ddj Xm as well as ud(j_i) m+ i has non-zero components and the others have the value of 
zero. This linear simultaneous equation is solved by conjugate gradient squared method (Dongarra 
et al. 1991). If we write equation( A.2) as A%jj = p, the algorithm is as follows: 



1. Prepare the initial guess ipo 
f = p- Atp Q 

Po = e = f 
I = 

2. while ||rj|| > e||p|| do 

(a) ai = (r ,fi)/(f ,Api) 



(b) h 



i+i 



ei - aiApi 
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(c) h+i =h~ aiA(ei + h l+1 ) 

(d) ip t+1 =xi + ai{ei + h t+1 ) 

(e) A = (r ,n + i)/(f ,n) 

(f) e m = ri+i + A^+i 

(g) Pl+i = ei+i + P{h+i + PiPi) 

(h) 1 = 1 + 1. 

The convergence of this method is improved much by preconditioning. In the present paper, we apply 
the incomplete LU-decomposition. The incomplete LU-decomposition is a method that matrix A is 
decomposed as LDU — R, where L, D, U, and R are a lower triangle matrix, a diagonal matrix, an 
upper triangle matrix, and a residual matrix. Instead of Aif) = p, we apply the CGS method to the 
reduced equation of 



{LDU)- l AiP = (LDU)- l p 



(A.3) 



Since (LDU)~ 1 A is nearly equal to the identity matrix / and the eigen values of this matrix gather 
around unity, the convergence of the CGS method is much accelerated. The way to decompose 
A = LDU — R is arbitrary. Here, we use the method by Meijerink & van der Vorst (1977) and 
Gustafsson (1978). Describing briefly, L, D, and U are chosen as 



( iM 

db 2 



\ 







dd m 
dc m+ i dd m+ i 

dc m+ 2 dd m+ 2 



db m+ i l/d m+ i 



V 



(A.4) 



db mn l/d mn J 



( k 







D 



(A.5) 







dmn—1 



dmn / 
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^ 1/di ubi udi 
l/d 2 ub 2 



U 



uc± 
ud 2 



UC 2 



1/d, 



m+1 



uh, 



m+1 







v 



(A.6) 



If we choose the diagonal element of D as 

dT 1 = a k - db k d k _iub k -i - dc k d 



k—m'U'C-k—m 



dd k d 



k a k-m+l u d k - m +l j 



(A.7) 



the diagonal part of the matrix A becomes equal to that of LDU . However, when the matrix LDU 
is operated on it generates extra terms (fill-in) such as 

db k d k ^iuc k -iip k 

+m— 1 

+ dc k d k _ m ub k ^ m ip k 

—m+1 

(A.8) 

as well as the ordinary terms: 

a k ^k + ub k ip k+1 + dbk^k-! + uc k tpk+m + dc k ipk-m + ud k 4>k+m-i + dd k ip k - m+1 . (A.9) 
To suppress the effect of these extra terms, the diagonal term d is modified as 

d' 1 = a k - db k d k _iub k _i- dc k d k _ m uc k - m - dd k d k _ m+l ud k _ m+ i-^{db k d k ^ 1 uc k _ 1 +dc k d k ^ m ub k _ m ), 

(A.10) 

where the factor 7 is chosen as the way by Gustafsson (1978). Using the LDU-decomposition, in the 
CGS algorithm the product of Av is replaced with (LDU)~ 1 Av. The operation multiplying (LDU) -1 
is done as follows: y = (LDU)~ 1 x is equivalent to x = (LDU)y, i.e., 



Lz = x 
DUy = z. 

The equation( A.ll| ) is solved as: for k = 1, 2, ■ • • , mn 

Zk = d k (x k - db k z k _i - dc k x k - m - dd k x k ^ m+ i), 
and then equation ( |A.12 ) is solved as: for k = mn, mn — 1, • • • , 1 

Vk = z k - d k (ub k y k+ i + uc k y k+m + ud k y k+m -i). 



(A.ll) 
(A.12) 

(A.13) 
(A.14) 



Necessary CPU time is proportional to the cell number, mn, and to loop count. A numerical 
experiment was done and we compared CGS with MILUCGS. The adopted test problem is as follows: 
the gravitational potential for a uniform cylinder is solved with square grids of m = n = N. The same 
problem was solved by a supercomputer Hitachi S-820/80, a general purpose main frame machine 
Hitachi M-682H, and a workstation Sun Sparc IPX. The allowable maximum relative error is chosen 
as e = 1CP 7 and calculations are done using 64 bits. As a result, about 80 x (A^/200) cycles are 
necessary to achieve a solution for MILUCGS, while the simple CGS requires 410 x (A^/200). Since 
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the number of cycles necessary for convergence is proportional to N, CPU time is proportional to iV 3 
irrespective with the preconditioning or not. 

In table Al, we summarize the CPU time necessary for solving the Poisson equation by respective 
machines. Comparing CPU time, the preconditioning reduces the CPU time about ^1/3 for scalar 
machines (M-682 and Sparc IPX). Since the amount of operations in MILUCGS is larger than that 
of CGS, the factor is not equal to ~ 1/5 (the ratio of necessary loops) but to <l/3. However, this 
clearly shows the advantage of MILUCGS over the simple CGS. The difference for the supercomputer 



is relatively small, for in order to vectorize the operation of multiplying (LDU) [equations (A. 11 



A. 12 )1 a more complicated algorithm with list vectors is necessary. 
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Figure Captions 

Fig.l: The initial radial density distributions. A solid line corresponds to models A and B, i.e., a 
magnetized cloud with a = 1. A dashed line shows the density of model C which has no magnetic 
fields a = 0. Both correspond to the models with the center-to-surface density ratio, F = 100. 
The dotted part of the curves represents the ambient medium / < 1 which is assumed to have no 
effect as a source of gravity. 

Fig. 2: The time evolutions of density (left panel) and magnetic fields (right panel) for model A. The 
z-axis and r-axis are directed horizontally and vertically, respectively, in contrast to an ordinary 
fashion. The number attached to the density contour lines represent the logarithm of the density: 
log 10 p. The step of the contour lines is taken constant = 0.25. Four snapshots, i.e., t = 2.4 (a), 
t = 3.6 (b), t = 4.9 (c), and t = 5.6 (d) are shown. In t = 3.6, no prominent fragmentations 
appear. In the stage shown in (c), a high-density region elongating in the z-direction appears. 
After (c), contraction along the symmetry axis proceeds. The final structure is a disk contracting 
towards the center seen in (d). The disk is perpendicular to the symmetric axis of the cylinder. 

Fig. 3: The structure of the cloud for model B. The upper-left panel shows the density distribution 
(a) and the upper-right does the magnetic field lines (b). It is shown that a fragment contracts 
and forms a disk perpendicular to the symmetric axis. Cross-cut views along the z-axis (c) and 
the r-axis (d) are plotted. In panel c, a logarithmic plot of p(z, 0) and B z (z, 0) and a linear plot of 
v z (z,0)/W are illustrated. Similarly, in panel d, a logarithmic plot of p(z c ,r) and B z (z c ,r) and a 
linear plot of v r (z c , r) /10 are shown, where z c represents the position of the center of the fragment 
= l z /2. This snapshot corresponds to the state of t = 1.26. 

Fig. 4: Time evolutions of the relative density enhancement, e) max , for models A and B. This is defined 
as <5 max = p(r = 0)| max / < p(r = 0) > —1. The curve for model B is similar to the late phase (t>4) 
of the model A. This shows that the structure and evolution of the fragment can be calculated 
from the initial state which has the most unstable perturbation with a finite amplitude. 

Fig. 5: The time evolution of density for Model C. In this model, the fragmentation in non-magnetic 
cloud is studied. Since the magnetic field, which has the effect of supporting cloud laterally, is not 
included, the radial extent of the cloud is thinner than that of models A and B. Two snapshots of 
density distribution at t = 1.2 (a) and t = 1.4 (b) are plotted. 

Fig. 6: The cross-cut views of the cloud of models B (a) and C (b). The left panel shows the 
variation of p(z, 0) against log \z — z c \, and the right panel shows that of p(z c , r) against logr. Two 
characteristic solutions [ Q3.4|) , and ( |3.6| )1 are also shown in right panels. Five snapshots at t = 0.3, 
0.6, 0.9, 1.2, 1.4 are plotted. 

Fig.7: The density distributions of models DB (a), DD (b), and DC (c). Model DB (a = 1) has the 
same parameter as model B except for the density contrast between the center and the surface, 
F = 10. Model DD corresponds to the cloud with stronger magnetic fields (a = 4). Model DC 
(a = 0) has the same parameter as model C except for the density contrast between the center 
and the surface, F = 10. For panel c, comparing how the density contour lines are running with 
that for Fig. 5, it can be seen that the density distributions are very similar with each other. Since 

— 1/2 

the length of numerical box in the z-direction l z is taken as proportional to p c , it is seen that 
the fragments of models C and F have an almost similar structure, if the size is scaled with the 
"scale-height" at the center, c s /(47rG/) c ) 1//2 instead of H. The number attached to the contour 
line, k, indicates that the value of the contour is equal to p = 10 fc / 4 . 
Fig. 8: The structure of the fragment of model CB2 at t = 0.705. Physical parameters of this model 
are chosen identical with model B but this model has 4.28 times higher spatial resolution than 
model B. In this model a perturbation with a larger amplitude <5j n i = 10 _1 is assumed initially 
to reduce a computing time. This time corresponds to 1.26 for model CB with <5j n i = 10~ 2 . The 
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panel a shows the density and velocity fields and b shows the magnetic field lines. To see the 
structure near the fragment clearly, in panel a, only a region of < z < l z /2 and r < l z /2 is 
plotted. Velocity vectors are plotted every 16 grids. Panels c and d are respectively, the cross-cut 
views along the z-axis and r-axis. A logarithmic plot of p(z, 0) and B z (z, 0) and a linear plot of 
v z (z, 0) are illustrated against log\z — z c \ in panel c. Similarly, in panel d, a logarithmic plot of 
p(z c ,r) and B z (z c ,r) and a linear plot of v r (z c ,r) are shown against logr. 
Fig. 9: The structure of the fragment formed in a cylindrical cloud threaded by a uniform magnetic 
field (model UB). The initial radial distribution of the density is as the same as model C (non- 
magnetic model), since the magnetic fields do not play a role in supporting the cloud. The density 
contour (a) and the magnetic fields (b) are illustrated. Panels c and d are, respectively, the cross- 
cut views along the z-axis and the r-axis. That is, a logarithmic plot of p(z, 0) and B z (z, 0) and a 
linear plot of v z (z, 0) are illustrated in panel c. Similarly, in panel d, a logarithmic plot of p(z c , r) 
and B z (z c ,r) and a linear plot of v r (z c ,r) are shown. This snapshot corresponds to the state of 
t = 1.26. 



